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ABSTRACT 

We model the surface magnetic field and open flux of the Sun from 1913 to 
1986 using a surface flux transport model, which includes the observed cycle- 
to-cycle variation of sunspot group tilts. The model reproduces the empirically 
derived time evolution of the solar open magnetic flux, and the reversal times of 
the polar fields. We find that both the polar field and the axial dipole moment 
resulting from this model around cycle minimum correlate with the strength of 
the following cycle. 

Subject headings: 



1. Introduction 



The evolution of the large scale magnetic field at the surface of the Sun can be modeled 
using a two-dimensional surface flux tr ansport rnodel where the magnetic fields undergo a 
random walk due to supergr anular flows flLeightonlll9641). are advected by differ ential rotation 
and meridional circulation ( lBabcocklll96ll : lLeightonlll964j: ISheelev et al.lll983l). and are sub- 



ject t o a slow decay from three-dimensional processes (jSchrijver et al 



2002 



Baumann et al. 



20061 ) . The surface f ield can be extrap olated out into the hehosphere, including the re- 
gion near the Earth ( Wang et al.l 120001 ). This makes the historic record of the magnetic 
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environment of the Earth, as manifested in the geomagnetic perturbation indices, a valu- 
able constraint on the variation of the magnetic fields of the Sun. The use of the surface 



Mackav et al. 


2002; 


Wane et al. 


2002b: 



The flux transport model has been refined over time: early work assum ed time-independe nt 



flows and did not in clude the effects of radial diffusion. It was found (iLean et al.l 12002 



P 



Schrijver et al.ll2002l ) that the observed variation in the cycle amplitudes would then lead to 
a secular drift of the polar field, in contradiction to observations. Two ways of extending 
the model were considered to cope with this problem: (1) making the meridional velocity 



time dependent (IWang et al. 



2002af). and f2) assuming that the poloida l field decays with 



a timescale of about 5 years (jSchrijver et al.ll2002l : iBaumann et al.ll2006l ). Here we consider 



a third possibility, namely that the tilt angle of the groups, the subject of Joy's law, varies 
from cycle to cycle. 

Most surface fiux transport studies have used fixed di f ferentia l rot ation and mer i diona l 
fiows. Two of the exceptions are the studies of IWang et al.l (l2002a[) andlDikpati et al.l (l2004j ) 
who considered cycle-to-cycle changes of the meridional flow. IWang et al.l (l2002al ) suggested 
that a suitably varying time-dependent meridional flow would prevent a secular drift of the 
polar flelds and thus allows them to reverse every cycle despite large variations in the cycle 
amplitudes. 



Our study is motivated by the recent flnding of iDasi-Espuig et al.l (|2010[ ) that the tilt 
angles of sunspot groups from the Mount Wilson Observatory and Kodaikanal observations 
( iHoward et al. l llQsi l ll999l : ISivaraman et al. show a cycle-to-cycle variation of Joy's law. 



Further they showed that the average tilt angle is negatively correlated with the strength 
of the cycle, i.e., the tilt angle is smaller for stronger cycles. A reduced tilt angle entails 
a smaller latitudinal separation between opposite polarity spots within a group, leading to 
reduced advec tion and diffusion of followin g polarity magnetic flux towards the poles during 
strong cycles (ICameron &: Schiisslerl 120071 ). 



Here we include the observed tilt angle variations as input to the flux transport model. 
Doing so requires us to reconsider the various parameters which go into the model (within 
the range constrained by observation s). In addition, we t entatively consider the effect of the 
observed inflows into active regions jHaber et al.l 120041 : iHindman et~all bood: iKomm et al.l 



20071 ) which cause a reduced escape of flux from active regions (jPe Rosa fc Schrijverll2006l ). 



The paper is organized as follows: Section 2 describes the flux transport model, including 
a brief discussion of how well the model parameters are observationally known. Section 3 
outlines the observations to which the model's results are compared in order to further 
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constrain the parameters and test the model. In Section 4 we give the parameters for our 
reference case. A brief parameter study, concentrating on the quahtative changes which 
occur as the different parameters are varied, is presented in Section 5. Our conclusions are 
given in Section 6. 



2. Flux transport model 



The surface flux transport model describes the passive transport of the radial compo- 
nent of the magnetic field, B, on the solar surface under the effects of differential rotation. 



Q, meridional flow, v, and sur 



diffusion (iDeVore et al. 



Schrijver et al. 



2002 



1985 



ace diffusivity, r? h, whilst gradually decaying owing to radial 



Baumann et al. 



Sheeley Jr et al.lll985t IWang et al. 



1989 



Mackav et all 12000 



20oX A source term, S'(A,0, t), describes the emer- 
gence of new flux as a function of latitude. A, and longitude, 0. The governing equation 
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where Z) is a linear operator describing t he decay due to radial d iffusion with radial diffusivity 
rjr. For D we adopt the form used in iBauraann et al.l ( 2006 ). We use the time averaged 
(synodic) differential rotation profile given by lSnodgrass (1983): f2(A) = 13.38 — 2.30 sin^ A — 
1.62 sin'^ A (in °/day). For t he time-averaged meridional flow we use the same proflle as 
van Ballegooijen et al.l ( 1l998l ). i.e. 



11 sin(2.4A) ms 




-1 



where|A| < 75° 
otherwise. 



(2) 



The two remaining parameters of the flux transport model are the horizontal and radial 
diffusivities, rjn and rjr- The r esults of various at t empts to measure rin from observation 
are summarized in Table 6.2 of ISchrijver fc Zwaaru (120001 ). The values obtained from cross- 
correlation and object-tracking methods fall in the range 100 — 300 km^s"^. We have used 
7]H = 250 km^s"^ for our reference value in Section 4; this value lies within the range of the 
observations, but we also consider the effect of varying it in Section 5. 
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Much less is known about rjr- This term was introduced by iBaumann et al.l (120061 ) to 
account for the 3D radial diffusion of the magnetic field and to obtain regularly reversing 
polar fields for cycles of varying amplitude in the absence of variations of the meridional 
flow. Its physical motivation is that the Sun's magnetic field is three-dimensional and thus 
has more modes of decay than are captured by the two-dimensional surface diffusion. We 
find here that the results with rjr = match the observations well (including having the 
polar fields reverse each cycle) when we include the observed tilt angle variations. We thus 
take ?7r = as our reference value and consider other values in Section 5. 



For the source t erm S (X,(l),t) in Equation [T] we follow Ivan Ballegooijen et al.l (Il998l ) 



and iBaumann et al.l (|2004l ) and consider new flux to emerge in the form of of opposite 
polarity patches. The positive-polarity patch is centered on latitude A+ and longitude 0+, 
the negative patch at (A_, 0_). The field of each new bipole is given hj B = — B^ with 

5±(A,0)=5^ax(^^) exp(2[l-cos(/3±(A,0))/52]), (3) 

where /3±(A,(/)) are the heliocentric angles between and (A±,(/)±), respectively and 

A/5 = /3+(A_,0-_) is the separation between the two polarities, and 5 = 4° is the size of 
the individual polarity patches. For the purposes of comparing the flux transport simula- 
tions with observations it is necessary to connect S closely to the actual observations. We 
use sunspot group areas and locations corresponding to their time of maximum area from 
http://solarscience.msf c.nasa.gov/greenwch.shtmil (based on the Greenwich photo- 



heliographic maps from 1874 to 1976 and USAF/NOAA SOON data thereafter) as proxies 
for emerging flux. 

The Greenwich/USAF/NOAA record contains the locations and areas of sunspots groups, 
but no magnetic polarity information. We use the location and areas to construct bipo- 
lar magnetic regions with the form described by Equation 3. The location of the bipoles, 
(A±, 0±), in the northern hemisphere are given by 

A± = A„± (-l)"0.5A/3sina (4) 
0± = (-l)"0.5A/3cosa(cosA)"\ (5) 

and those in the southern hemisphere by 

A± = A™± (-l)"0.5A/3sina (6) 
0± = (-l)"0.5A/3cosa(cosA)"\ (7) 

Here (A^, <Pm) is the central location of the group from the Greenwich/USAF/NOAA record, 
a the tilt angle with respect to the azimuthal direction, and n the cycle number. The 
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1/2 

separation, A/3, between the two polarities is taken to be A/3 = 0.45^4^ where is the 
total area of the active region. We estimate the total flux of an active region by considering 
its total area, Aji, to be the sum of the area covered by the sunspots Ag and the plage Ap 
using the observed relationship 



A 



R 



As + A„ = A, + 414 + 21 A - 0.0036^^ 



(IChapman et al.l 119971 ). where all areas are measured in millionths of a solar hemisphere. 
The coefficient 0.45 was determined by using the sunspot group data from Mount Wilson 
(covering the period from 1917 to 19 85) and Kodai k anal (cove ring the period fr o m 19 06 to 
1987). The data sets are described in lHoward et al.l (119841 ) and ISivaraman et al.l (119991 ) and 
are available from http://ngdc.noaa.gov/stp/SOLAR/ftpsunspotregions.html. Both data 
sets include umbral areas Au and separations A/3 between the "centers of mass" of the 
leading a nd following spots. We converted the umbral areas to sunspot areas using the 
results of iBrandt et al.l (1l990[ l. and from there to Ar using Equation 8. Figure [1] shows the 
average (over 7-degree bins) of A/3 from each data set together with the fit curve. 

The next step is to specify the tilt angle, a, including its cycle-to-cycle variations. As 
noted in iBaumann et al.l (120041 ) the polar fields are essentially proportional to a, so that 
these variations might strongly affect the results. We use the tilt angles provided by the 
sunspot group data from Mount Wilson and from Kodaikanal. Since these observations 
cover only part of the sunspot groups in the combined Greenwich/USA/NOAA dataset that 
we use for the source, we take cycle-averaged properties for the tilt angle as a function of 
latitude. The asterisks in Figure [2] show the binned cycle averages of the tilt angle weighted 
with the group areas. We fit the data from each cycle to the form a = Tn\/JX\, where n is 
the cycle number. We calculate T„ by 

T„ = ^^^-^ (9) 



where the summation is over all spot groups in cycle n; Ai is area of the i-th spot group and Aj 
is its latitude. If the tilt angle of each group of cycle n is written as ctj = Tn\/JXi\+ei, then the 
above estimate for T„ implies Yli ^i^i = i-^-, the area- weighted sum of the deviations from 
the fit curve is zero. We see in Figure [2] that the square-root profile produces a reasonable 
fit. Furthermore cycle 19 is systematically low across a broad range of latitudes. Figure [3] 
shows T„ for both the Mount Wilson and Kodaikanal data sets. 



Observed local i zed in flows associated with active regions (IHaber et al.ll2004j : lHindman et al 



2004 



Komm et al.l 120071) a re important for the evolution of the surface magnetic fields 
(iDe Rosa &: Schrijverl 120061 ) and should be included in the model. The effect of the in- 



flows is to reduce the rate of expansion of the active region flux. This reduces the latitudinal 
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separation of the polarities and thus the amount of net fl ux which can migrate t o the poles. 
Additionally the large-sc a le inflo ws into the activity belt (iGizon &: Rempelll2008l ) affects the 
polar fields. IJiang et al.l ( 12 01 Obi ) have quantitatively studied the effects of such inflows on 
the polar fields, finding that a 5 m/s global inflow such as was observed during cycle 23 re- 
duces the polar fields by 18%. The flows into indi vidual active regions are str onger than this 
global-scale inflow and thus have a bigger effect (jPe Rosa fc Schrijverl 120061 ). In this paper 
we have tentatively introduced a time-independent parameter, g, which we use to scale the 
observed tilt angles of the sunspot groups. The main effect of g is to scale the amount of 
flux which reaches the pole against the amount of flux which emerges in the sunspot groups. 
We comment that the introduction of this parameter also deals with the uncertainty of how 
the observed group tilt angles (which are based upon the position of sunspots in continuum 
intensity images) are related to the tilt angles between the opposite polarity patches of the 
active region. 

Finally we determine -Bmax in Equation 3 by matching the unsigned total observed flux 
from the Mount Wilson and Wilcox observatories with the simulation results. 

For the initial niagneti c field distribution at the start of the simulation we follow 
van Ballegooijen et al.l (|l998l ) and use an axisymmetric solution to Equation 1 with B = ±Bq 
at the poles, and which evolves almost entirely on the slow diffusive timescale: 



B 



sign(A)-Bo exp 
sign(A)5o 



-llms ^xR, 



2AriH 



1 + COS2.4A) 



if |A| < 75° 
otherwise. 



(10) 



In the absence of sources, the evolution of this initial condition is dominated by the slow 
decay of the global field (with an e- folding time of approximately 4000 years). This choice 
for the functional form of the initial condition is arbitrary. The time between flux emerging 
and its reaching the poles is on the order of 5 years, so the polar fields for the first 5 years 
or so are strongly affected by the chosen form for the initial condition. We therefore exclude 
the first polar field maximum from our analysis of the results. 



3. Open flux and timing of the polar reversals 

The model described above requires emerging bipolar magnetic region areas, locations, 
and tilt angles as input. We have these data for the period between 1913 and 1986. The 
flux transport model gives as output the radial component of the magnetic field on the 
solar surface. Throughout most of the period covered by the simulations, observational 
magnetogram data are unavailable for comparing against the results of the simulations. We 
therefore consider the Sun's open flux, -Fopen, which has been inferred from the aa-index of 
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Fig. 2. — Area weighted tilt angle Aiai)/{J2i ^i) a function of sunspot group latitude 
for cycles 15 through 21 on the basis of the combined Mt. Wilson and Kodaikanal datasets. 
The asterisks represent the average values for 5 degree bins for each cycle. The curves 
represent fits of the form a = Tn^^/\X\ for each cycle n. The blue asterisks and curve indicate 
cycle 19. 
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geomagnetic var iations (and its extensions) from 1842 onwards (ILockwood &: Stamperlll999 
LockwoodI I2OO3I ) ■ To obtain Fopen from the simulat ion, we take the surface d is tribution of 
B an d use the current sheet source surface model ( Zhao &: Hoeksema 1995a lb: Zhao et al. 
2002[) to extrapolate the solar surface field out into the heliosphere (see Schiissler fc Baumann 
20061 : I Jiang et al.ll2010al ). The results of the extrapolation are dependent on the assumed 
value of the 'cusp radius', -Rcusp, which is the radial distance beyond which all the field lines 
are open. 

We also compare the sir nulation results agains t the timing of the polar field reversals, 
which have been inferred by iMakarov et al.l (120031 ) from polar filament observations from 
1870 to 2001. 



Reference case 



We here give the results for the parameter set rj^ = 250 km^s~^, g = 0.7, Bq = —10.2 
G, rjr = and -Rcusp = 1.55i?Q. The value -Bmax = 374 G was found by matching to the total 
observed unsigned magnetic flux from the Mount Wilson and Wilcox Solar Observatories (see 
Figure [H) . With these parameters, the model reproduces well the open flux -Fopen inferred by 
LockwoodI (120031 ) as shown in Figure [5l This applies to the phases as well as the amplitudes 
of both the maxima and minima of the inferred open flux. The corresponding evolution of 
the polar field (defined as the average field above ±75° latitude) and axial dipole moment 
are shown in the upper panel of Figure [61 The polar field closely follows the axial dipole 
moment, with a delay of several years. This is understandable as the dipole moment reacts 
more quickly to flux transport across the equator, which then takes several additional years 
to reach the polar latitudes (> 75°). The simulated polar fields reverse for all cycles. Without 
the cycle-dependent variations of the tilt angle the weak cycle 20 would ha ve been unable 



to offs et the polar field after cycle 19. It was this type of problem which led lSchrijver et al. 
(I2OO2I ) and iBaumann et al.l (120061 ) to introduce a decay term. Here we achieve a good 



agreement with the observations even without su ch a term because of the anti-correlation 
in the observed tilt angles and cycle strengths (jPasi-Espuig et al.l I2OIOI ). The asterisks 
in the upper panel o f Figure 6 indicate the timings of the polar reversals as derived by 
Makarov et al.l ( 20031 ) from Ha polar filament data. The reversal times are reasonably well 
reproduced, except for the first reversal which is still affected by arbitrary form of the initial 
condition. 

The maximum of the dipole moment during the activity minimum between cycles 19 
and 20 is much lower than that between cycles 18 and 19. This is primarily because the 
average tilt angle of cycle 19 was substantially lower than that of cycle 18 (see Figures 2 and 
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3). On the other hand, the dipole moment between cycles 20 and 21 is again high, because 
the average the tilt angles of cycle 20 was high. The lower panel of Figure [6] presents the time 
evolution of the unsigned north-south averaged polar field from the simulation as well as the 
observed strength of the solar activity. Ignoring the first cycle, which is strongly affected 
by the initial condition, the polar field can be seen to anticipate the peaks of solar activity. 
We have quantified this relationship by calculating the correlation between the strength 
of peaks the (unsigned north-south averaged) polar field and the strength of the adjacent 
activity maxima. As can be seen in Figure [TJ the maxima of the polar field between cycles 
n and n + 1 reflect the amplitude of cycle n + 1 much better than cycle n. The correlation 
coefficient for the relationship between the polar field and the ampli tude of the next cycle is 
0.85. This is in contrast to the result of ICameron &: Schiisslerl (120071 ) who did not consider a 
cycle-to-cycle dependent tilt angle; they found that the polar fields and strength of the polar 
field closely follow the previous cycle. This indicates the importance of the cycle-to-cycle 
variations in the tilt angle. 

Some observational evidence concerning the correlation between the polar field and 
th e strength of the prev i ous and subsequent cycles has been previo u sly consider e d, e.g , 
bv lSchatten et all dlOTs I): llavden et all fll99lh : ISvalgaard et~aD J2005[ l: Ijiang et all feoOTh . 
The existence of a correlation between the polar field and the strength of the next cycle is 
evidence in favour of a Bacock-Leighton-type dynamo. Within the co ntext of such dynam os, 
the correlation constrains the subsurface dynamics (see, for example, lYeates et al.ll2008l ). 



5. Parameter dependence 

In the previous section we showed that a good representation of the empirically deter- 
mined open flux -Fopen could be achieved with the parameter values of the reference model. 
Figure [8] shows the effect of separately varying the initial field strength, Bq, the surface 
diffusivity, rju, radial diffusivity, t]^, and cusp surface height, -Rcusp- In all panels except the 
lower right we have kept -Bmax constant: in this panel we recalibrated -Bmax to account for 
the change in the total unsigned fiux resulting from the change in rj^. 

The upper left panel of Figure |8] shows the effect of changing Bq. Since Bq describes 
the initial polar field, varying this value leads to an offset in the strength of the axial 
dipole moment, which persists throughout the simulation when 77^ = 0. This offset results 
in alternating cycles having either stronger or weaker axial dipole moments depending on 
whether or not they have the same sign of dipole moment as that of the initial state. Near 
activity minima it is the lower order axial moments which dominate -Fopen so that the minima 
alternately become higher and lower. 
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The upper right panel of Figure [8] shows the effect of increasing rjH'- the minima of 
-Fopcn are shifted upwards, while the maxima are not substantially affected. The explanation 
for the upward shift is that tih determines the amount of flux which crosses the equator and 
thus directly influences the axial dipole moment. There is also a weak but noticeable 22-year 
component, with the minima of alternating cycles being weaker. This 22- year component is 
present because we have not recalibrated Bq. 

The middle left panel of Figure [S] shows the effect of increasing rj^. The enhanced decay 
of the field not only reduces -Fopen around the minima but also during the rise phase a cycle. 
This leads to too low minima and a delay of the rising phase. There is again a 22-year 
component because Bq has not been recalibrated. 

The middle right panel shows the effect of varying the tilt angle reduction factor, g, from 
0.7 to 1, which modifies the magnitude of the polar fields and axial dipole moment. The 
signature is therefore an increase in the magnitude of the changes in the dipole moment (and 
thus -Fopen) its minima, so that the effect almost cancels after two cycles. This also produces 
a strong 22 year periodicity in the minima. We comment that g = 0.7 is required to obtain 
the correct ratio between the maxima and minima of the open flux, as it essentially scales 
the low-order axial multipoles whilst barely affecting the equatorial multipoles. Introducing 
g does not affect whether or not the polar fields reverse - the 22 year periodicity, when g is 
varied in isolation, can be removed by an appropriate choice of Bq. 

The lower left panel shows that increasing -Rcusp in isolation weakens -Fopen- The effect is 
strongest during the maxima as it preferentially reduces the contribution from higher order 
multipoles. The influence is thus qualitatively different from that of the other parameters in 
that it changes the relative contributions of the different multipoles. 

In the panels discussed so far we have kept -Bmax, the scaling factor for the total flux of 
newly emerging BMRs, constant. Varying rjj. as was done in the middle left panel changes 
the total amount of unsigned flux, and so affects the calibration of -Bmax- In the bottom 
right panel we therefore show the effect of a change in 77^ together with the corresponding 
change in -Bmax- Since the entire system is linear in -Bmax, changing -Bmax nierely rescales the 
result - hence the result in the lower right panel is just a scaled version of the result shown 
middle left panel. We note that varying rjH also affects the calibration. 

This brief study of the effect of varying the parameters illustrates the kind of changes 
which occur. However, it does not rule out other choices for the parameters which also could 
provide a good fit to the observations. In particular we do not claim that non-zero values 
of rjr are excluded, although we can say that, at least for cycles 15-21, a good fit to the 
observations does not require a non-zero rjr. 



6. Conclusions 



The surface flux transport model including the effect of a cycle- dependent variation of 
the tilt angles of sunspot groups reproduces the major features of the observationally inferred 
open flux and the timing of the polar field reversals from 1913 to 1986 (the period for which 
we have the tilt angle data) . The reversal of polar fields after strong cycles can be explained 
by the observed anti-correlation between the active region tilt angle and cycle amplitude, so 
that no additional decay by radial diffusion was required to achieve this result. 

When the observed tilt angle are used, the polar field maxima from the model are 
correlated with the strength of the following cycle. This correlation is hkely to be present 
independent of the parameter choices, provided the model reproduces the minima of the 
inferred open fiux. The correlation suggests that the polar fields are an important ingredient 
of the solar dynamo process, which is consistent with Babcock-Leighton-type models. The 
cycle-to-cycle variation of Joy's law might play a role in the nonlinear modulation of the 
solar dynamo. 
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Fig. 3. — The cycle dependence of the factor T„ in the fit curves of Figure [2l Shown are the 
values derived from the Mt. Wilson data (blue) the Kodaikanal data (red). 
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Fig. 4. — Average unsigned magnetic field from tfie simulation with -Bmax = 374 G (solid 
curve) compared with observations from the Mount Wilson (plus signs) and Wilcox Solar 
Observatories (squares). 
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Fig. 5. — Open flux inferred from tfie aa-index (ILockwoodl l2003l black) and result from 
the surface flux transport model and extrapolation (red) for the reference case. The dashed 
vertical black line refers to the time beyond which we have no observed values of the tilt 
angle. 
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Fig. 6. — Upper panel: axial dipole moment (black) together with the north (solid red) 
and south (dashed red) polar field strengths from the fiux tran sport simulation . The black 
asterisks indicate the inferred times of polar field reversals from iMakarov et al.l (120031 ) . The 
dashed vertical black line refers to the time beyond which we have no observed values of the 
tilt angle. Lower panel: average of the unsigned polar field strength from the fiux transport 
model (red) and observed sunspot (black) . 
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Fig. 7. — Correlation between the observed maximum polar field from the flux transport 
model and and the preceding (upper panel) and subsequent (lower panel) cycle's maximum 
sunspot area. The asterisks represent the values for individual cycles: the solid line is a 
linear fit. The maximum sunspot area is uncorrelated (r = —0.1) to the subsequent polar 
field maximum, but strongly correlated (r = 0.85) to the preceding polar field maximum. 
The first maximum of the polar field has been omitted from the analysis since it is affected 
by the arbitrariness of the initial field. 
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Fig. 8. — Effect on -Fopen of varying the parameters of the flux transport model. The black 
curve in each case is the observationally inferred open flux while the red curve shows the 
value from the simulation. For all other parameters the reference values were used. The 
parameters varied were: 

upper left: initial field strength Bq (—10.2 G — )■ —8.5 G), 

right: surface diffusivity t]h (250 km^s^^ — > 500 km^s^^), 
middle left: radial diffusivity rjr (0 — ?■ 100 km^s~^ ), 

right: tilt angle reduction due to inflows g (0.7— t-I) 
lower left: cusp surface radius -Rcusp (1.55-Rq — 1.8-Rq) 

right: radial diffusivity r]r (0 ^ 100 km^s'^) and B^^^{374 G 462 G) 



